Method for wavefield separation for dual-sensor data using kirchhoff-type datuming and migration

ABSTRACT

A first weighted integral operator is applied to dual-sensor data to extrapolate the dual-sensor data to a first position above an acquisition surface, generating extrapolated data. A second weighted integral operator is applied to the extrapolated data to extrapolate the extrapolated data to a second position, generating wavefield separated data. One of the integral operators is applied to a scaled combination of the dual sensor data.

CROSS-REFERENCES TO RELATED APPLICATIONS

Not Applicable

FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

Not Applicable

SEQUENCE LISTING, TABLE, OR COMPUTER LISTING

Not Applicable

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates generally to the field of geophysicalprospecting. More particularly, the invention relates to the field ofwavefield separation in dual-sensor data.

2. Description of the Related Art

In the oil and gas industry, geophysical prospecting is commonly used toaid in the search for and evaluation of subterranean formations.Geophysical prospecting techniques yield knowledge of the subsurfacestructure of the earth, which is useful for finding and extractingvaluable mineral resources, particularly hydrocarbon deposits such asoil and natural gas. A well-known technique of geophysical prospectingis a seismic survey. In a land-based seismic survey, a seismic signal isgenerated on or near the earth's surface and then travels downward intothe subsurface of the earth. In a marine seismic survey, the seismicsignal may also travel downward through a body of water overlying thesubsurface of the earth. Seismic energy sources are used to generate theseismic signal which, after propagating into the earth, is at leastpartially reflected by subsurface seismic reflectors. Such seismicreflectors typically are interfaces between subterranean formationshaving different elastic properties, specifically sound wave velocityand rock density, which lead to differences in acoustic impedance at theinterfaces. The reflected seismic energy is detected by seismic sensors(also called seismic receivers) at or near the surface of the earth, inan overlying body of water, or at known depths in boreholes andrecorded.

The resulting seismic data obtained in performing a seismic survey isprocessed to yield information relating to the geologic structure andproperties of the subterranean formations in the area being surveyed.The processed seismic data is processed for display and analysis ofpotential hydrocarbon content of these subterranean formations. The goalof seismic data processing is to extract from the seismic data as muchinformation as possible regarding the subterranean formations in orderto adequately image the geologic subsurface. In order to identifylocations in the Earth's subsurface where there is a probability forfinding petroleum accumulations, large sums of money are expended ingathering, processing, and interpreting seismic data. The process ofconstructing the reflector surfaces defining the subterranean earthlayers of interest from the recorded seismic data provides an image ofthe earth in depth or time.

The image of the structure of the Earth's subsurface is produced inorder to enable an interpreter to select locations with the greatestprobability of having petroleum accumulations. To verify the presence ofpetroleum, a well must be drilled. Drilling wells to determine whetherpetroleum deposits are present or not, is an extremely expensive andtime-consuming undertaking. For that reason, there is a continuing needto improve the processing and display of the seismic data, so as toproduce an image of the structure of the Earth's subsurface that willimprove the ability of an interpreter, whether the interpretation ismade by a computer or a human, to assess the probability that anaccumulation of petroleum exists at a particular location in the Earth'ssubsurface.

The appropriate seismic sources for generating the seismic signal inland seismic surveys may include explosives or vibrators. Marine seismicsurveys typically employ a submerged seismic source towed by a ship andperiodically activated to generate an acoustic wavefield. The seismicsource generating the wavefield may be of several types, including asmall explosive charge, an electric spark or arc, a marine vibrator,and, typically, a gun. The seismic source gun may be a water gun, avapor gun, and, most typically, an air gun. Typically, a marine seismicsource consists not of a single source element, but of aspatially-distributed array of source elements. This arrangement isparticularly true for air guns, currently the most common form of marineseismic source. In an air gun array, each air gun typically stores andquickly releases a different volume of highly compressed air, forming ashort-duration impulse.

The appropriate types of seismic sensors typically include particlevelocity sensors, particularly in land surveys, and water pressuresensors, particularly in marine surveys. Sometimes particle displacementsensors, particle acceleration sensors, or pressure gradient sensors areused in place of or in addition to particle velocity sensors. Particlevelocity sensors and water pressure sensors are commonly known in theart as geophones and hydrophones, respectively. Seismic sensors may bedeployed by themselves, but are more commonly deployed in sensor arrays.Additionally, pressure sensors and particle velocity sensors may bedeployed together in a marine survey, collocated in pairs or pairs ofspatial arrays.

A particle motion sensor, such as a geophone, has directionalsensitivity, whereas a pressure sensor, such as a hydrophone, does not.Accordingly, the upgoing wavefield signals detected by a geophone andhydrophone located close together will be in phase, while the downgoingwavefield signals will be recorded 180 degrees out of phase. Varioustechniques have been proposed for using this phase difference to reducethe spectral notches caused by the surface reflection and, if therecordings are made on the seafloor, to attenuate water borne multiples.It should be noted that an alternative to having the geophone andhydrophone co-located, is to have sufficient spatial density of sensorsso that the respective wavefields recorded by the hydrophone andgeophone can be interpolated or extrapolated to produce the twowavefields at the same location.

Conventional 3D marine seismic acquisition by towed streamer usuallyresults in asymmetrical spatial sampling and fold between inline andcross-line directions. The sampling density is denser in the inlinedirection (parallel to the towed streamers) than in the cross-linedirection (perpendicular to the towed streamers). The asymmetry is dueto a wider spacing between receivers in separate streamers than betweenreceivers in the same streamer. This asymmetry can lead to spatialaliasing of the sampling data in the cross-line direction. When thereceivers are deployed on the sea-bottom, the lateral distance betweenadjacent receivers is often larger than the typical receiver separationin a streamer. This means that spatial aliasing can be a problem fordata acquired with this type of receiver configuration as well. Thisinvention considers data acquired with hydrophones and geophones. Thegeophone can record one or more components of the particle velocityvector. Common are three-component geophones which record the completeparticle motion vector. Recording only the vertical component of theparticle velocity is also common. Hydrophone and geophone recordstogether are referred to as dual-sensor data in this invention.

Thus, a need exists for a method for wavefield separation fordual-sensor data that more efficiently corrects for spatial aliasing.Preferably, the method accomplishes this in the time-space domain whereirregular receiver separations can be taken into account and does notrequire seismic trace interpolation.

BRIEF SUMMARY OF THE INVENTION

The invention is a method for wavefield separation for dual-sensor data.A first weighted integral operator is applied to dual-sensor data toextrapolate the dual-sensor data to a first position above anacquisition surface, generating extrapolated data. A second weightedintegral operator is applied to the extrapolated data to extrapolate theextrapolated data to a second position, generating wavefield separateddata. One of the integral operators is applied to a scaled combinationof the dual sensor data.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention and its advantages may be more easily understood byreference to the following detailed description and the attacheddrawings, in which:

FIG. 1 is a schematic side view of a ray from the acquisition surface;

FIG. 2 is a section showing a traveltime operator on top of a part of acommon-shot pressure record;

FIG. 3 is a flowchart illustrating general embodiment of the inventionfor wavefield separation for dual-sensor data;

FIG. 4 is a flowchart illustrating a first embodiment of the inventionfor wavefield separation for dual-sensor data;

FIG. 5 is a flowchart illustrating a second embodiment of the inventionfor wavefield separation for dual-sensor data;

FIG. 6 is a schematic side view of a single ray corresponding to a timesample from a trace to be migrated;

FIG. 7 is a flowchart illustrating a third embodiment of the inventionfor wavefield separation for dual-sensor data;

While the invention will be described in connection with its preferredembodiments, it will be understood that the invention is not limited tothese. On the contrary, the invention is intended to cover allalternatives, modifications, and equivalents that may be included withinthe scope of the invention, as defined by the appended claims.

DETAILED DESCRIPTION OF THE INVENTION

The invention is a method for wavefield separation for dual-sensor data.In the following, three embodiments of the method of the invention forseparating recorded wavefields into their up-going and down-going partsare described. Wavefields are recorded with receivers disposed in aplurality of seismic streamers, also known as a 3D spread, or deployedon the sea-bottom. The receivers define the acquisition surface, whichis a horizontal plane in the ideal case of flat horizontal streamers ora flat horizontal sea-bottom. However, the acquisition surface need notbe flat, as this is not a limitation of the invention. At each receiverlocation on the acquisition surface, two wavefields are recorded orbuilt by means known in the art from more recorded wavefields. The twowavefields are the pressure field consisting of up-going and down-goingparts (hydrophone record) and the component of the particle velocityfield normal to the acquisition surface (geosensor record), againconsisting of up-going and down-going parts. The method of the inventionemploys a scaled combination of the two recorded wavefields in theintegral expressions that separate the wavefields.

Conventional methods consider only homogeneous media. In the method ofthe invention, the earth model can be laterally heterogeneous. However,the earth model can also be replaced by a smooth background model whichcorrectly represents the kinematics of wave propagation and is inaccordance to the high-frequency approximation. This approximation meansthat wave propagation can be described with high accuracy by zero orderray theory. All three embodiments exploit the principle of stationaryphase which automatically selects the correct obliquity factor to beapplied in the wavefield separation.

The down-going parts of the wavefield are due to the reflection of allup-going events at the sea-surface. In the following, a hypotheticalmodel where the sea-surface is absent will be considered. The down-goingenergy will be treated implicitly as up-going energy. This correspondsto a virtual reflection of the down-going energy at the acquisitionsurface without any amplitude changes, i.e. without changes to magnitudeand polarization. The only thing which happens mathematically is achange of sign of the vertical component of the slowness vectorassociated with a down-going wavefront. The space above the sea-surface,filled with air in reality, is assumed to be filled with water such thata smooth transition of all medium properties at the sea-surface isassured.

The first embodiment of the method of the invention is wavefieldseparation using Kirchhoff-type time domain or frequency domaindatuming. This first embodiment extrapolates a combination ofcommon-shot records of the recorded wavefields to a level above theacquisition surface. The combination of the recorded wavefields is suchthat they separate into up-going and down-going parts in this datumingstep. The separated wavefields are then back-propagated to the originallevel, generally to any level. The output of this first embodiment iscommon-shot ensembles of seismic traces, i.e. seismic pre-stack data ofseparated wavefields.

A pressure field p(x, ω) recorded at a point x on an acquisition surfaceS can be forward extrapolated in time to another point x^(R) above theacquisition surface using the following integral equation:

$\begin{matrix}{{{p\left( {x^{R},\omega} \right)} \cong {\frac{\omega}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}{p\left( {x,\omega} \right)}{\exp \left\lbrack {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right\rbrack}\ {S}}}}},} & (1)\end{matrix}$

where i is the imaginary unit and ω denotes angular frequency. Thetraveltime along the ray connecting x and x^(R) is given by τ(x, x^(R)).The term W(x, x^(R)) in Equation (1) is a weighting factor, which isgiven by

$\begin{matrix}{{{W\left( {x,x^{R}} \right)} = {\frac{1}{\rho (x)}{G_{0}\left( {x,x^{R}} \right)}\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{c(x)}}},} & (2)\end{matrix}$

where ρ(x) is the density and c(x) is the speed of sound in the medium.The factor G₀(x, x^(R)) in Equation (2) takes into account amplitudechanges due to geometrical spreading. FIG. 1 is a schematic side view ofa ray from the acquisition surface. θ(x, x^(R)) denotes the anglemeasured at x between the ray connecting points x and x^(R) and the unitnormal vector to the acquisition surface, n(x), which points towardsx^(R). The equation of motion relates the pressure field with thethree-component particle velocity field. The component of the particlevelocity field normal to the acquisition surface ν_(n)(x, ω) is derivedfrom a three-component particle velocity field measurement by means wellknown in the art. The relationship between pressure and particlevelocity fields allows using an integral equation analogous to Equation(1) for datuming the normal component of the particle velocity:

$\begin{matrix}{{v_{n}\left( {x^{R},\omega} \right)} \cong {\frac{\omega}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}{v_{n}\left( {x,\omega} \right)}{\exp \left\lbrack {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right\rbrack}\ {{S}.}}}}} & (3)\end{matrix}$

In a homogeneous medium with constant media properties ρ and c, an eventin the pressure record or in the normal component of the particlevelocity record can be separated into its up-going and down-goingcomponents using the following equations:

$\begin{matrix}{{{p^{up}\left( {x,\omega} \right)} = {\frac{1}{2}\left( {{p\left( {x,\omega} \right)} - {\frac{\rho \; c}{\cos \left\lbrack {\theta (x)} \right\rbrack}{v_{n}\left( {x,\omega} \right)}}} \right)}},} & (4) \\{{{p^{dn}\left( {x,\omega} \right)} = {\frac{1}{2}\left( {{p\left( {x,\omega} \right)} + {\frac{\rho \; c}{\cos \left\lbrack {\theta (x)} \right\rbrack}{v_{n}\left( {x,\omega} \right)}}} \right)}},} & (5) \\{{{v_{n}^{up}\left( {x,\omega} \right)} = {\frac{1}{2}\left( {{v_{n}\left( {x,\omega} \right)} - {\frac{\cos \left\lbrack {\theta (x)} \right\rbrack}{\rho \; c}{p\left( {x,\omega} \right)}}} \right)}},{and}} & (6) \\{{{v_{n}^{dn}\left( {x,\omega} \right)} = {\frac{1}{2}\left( {{v_{n}\left( {x,\omega} \right)} + {\frac{\cos \left\lbrack {\theta (x)} \right\rbrack}{\rho \; c}{p\left( {x,\omega} \right)}}} \right)}},} & (7)\end{matrix}$

respectively, as is well known in the art. Here, cos[θ(x)] denotes theangle between the propagation direction of the wavefront associated withthat event and the unit normal vector n, and ν_(n)(x, ω) is the recordedcomponent of the particle velocity field normal to the acquisitionsurface S. Equations (4) through (7) remain valid in slightlyinhomogeneous media within the validity range of the high-frequencyapproximation, in which a propagating wave behaves locally like a planewave, and, thus, the wavefield can be separated using the local mediumproperties ρ(x) and c(x). The angle θ(x) is initially unknown, but willbe automatically determined by application of the integral equation ofEquations (1) or (3), as the major contribution comes from a small zone(Fresnel zone) around the stationary point. At the stationary point, thetraveltime surface defined by τ(x, x^(R)) and the seismic event aretangent, i.e. both have the same dip, as is shown in FIG. 2.Consequently, when the correct local velocity is used,

cos[θ(x, x ^(R))]=cos[θ(x)].   (8)

The sum (integral) over contributions from outside a small range aroundthe stationary point is negligibly small. This localness of thecontributions remains valid when the wavefield is scaled by a smoothlyvarying function. Thus, only the contributions close to the stationarypoint and scaled by obliquity factors close to the correct obliquityfactor affect the results of the corresponding integral expressions.Thus, in a first embodiment of the method of the invention, theextrapolation of the up-going and down-going pressure fields and normalcomponents of particle velocity fields at point x^(R) are estimatedseparately by combining Equations (4) through (7) with Equation (1),yielding the following integral equations:

$\begin{matrix}{{{{p^{up}\left( {x^{R},\omega} \right)} \cong}\quad}{\quad {{\frac{\omega}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \left( {{p\left( {x,\omega} \right)} - {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,\omega} \right)}}} \right){\exp \left\lbrack {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right\rbrack}\ {S}}}},}}} & (9) \\{{{{p^{dn}\left( {x^{R},\omega} \right)} \cong}\quad}{\quad {{\frac{\omega}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \left( {{p\left( {x,\omega} \right)} + {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,\omega} \right)}}} \right){\exp \left\lbrack {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right\rbrack}\ {S}}}},}}} & (10) \\{{{{{v_{n}^{up}\left( {x^{R},\omega} \right)} \cong}\quad}{\quad {\frac{\omega}{2\pi}{\int_{S}{W\left( {x,x^{R}} \right)}}}\quad}\frac{1}{2} \left( {{v_{n}\left( {x,\omega} \right)} - {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,\omega} \right)}}} \right){\exp \left( {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right)}\ {S}},{and}} & (11) \\{{{{v_{n}^{dn}\left( {x^{R},\omega} \right)} \cong}\quad}{\quad {\frac{\omega}{2\pi}{\int_{S}{W\left( {x,x^{R}} \right)}}}\quad}\frac{1}{2} \left( {{v_{n}\left( {x,\omega} \right)} + {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,\omega} \right)}}} \right){\exp \left( {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right)}\ {{S}.}} & (12)\end{matrix}$

Using the well-known relationships between the frequency domain and thetime domain

$\begin{matrix}\left. {{\omega}\; {p\left( {x,\omega} \right)}}\leftrightarrow{\frac{\partial}{\partial t}{p\left( {x,t} \right)}\mspace{14mu} {and}} \right. & (13) \\{\left. {{p\left( {x,\omega} \right)}{\exp \left( {- {{\omega\tau}\left( {x,x^{R}} \right)}} \right)}}\leftrightarrow{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)} \right.,} & (14)\end{matrix}$

the time-domain extrapolation expressions corresponding to Equations (9)through (12) are obtained as:

$\begin{matrix}{{{p^{up}\left( {x^{R},t} \right)} \cong {\frac{1}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}\  - {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right){S}}}}},} & (15) \\{{{p^{dn}\left( {x^{R},t} \right)} \cong {\frac{1}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}\  + {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right){S}}}}},} & (16) \\{{{v_{n}^{up}\left( {x^{R},t} \right)} \cong {\frac{1}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}\  - {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right){S}}}}},{and}} & (17) \\{{{v_{n}^{dn}\left( {x^{R},t} \right)} \cong {\frac{1}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}\  + {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right){S}}}}},} & (18)\end{matrix}$

respectively. Either the time-domain extrapolation expressions inEquations (15) through (18) or the frequency-domain extrapolationexpressions in Equations (9) through (12) can be used for performingthis first datuming step in this first embodiment of the invention.

The maximum time dip in the output is determined by the chosen operatortaper in the first datuming step, i.e. by the maximum angle θ(x, x^(R))considered along the traveltime operator. The datuming distance thencontrols the number of traces in the aperture and the dip sampling asthe angle θ(x, x^(R)) generally decreases with increasing datumingdistance. FIG. 2 shows a traveltime operator on top of a part of apressure record. At the stationary point, the traveltime operator istangential to the seismic event. Both have the same dip. To provide theup-going and down-going fields at the original acquisition surface or alevel near by (e.g. for comparison with a different survey), a seconddatuming step is required, usually a backward propagation in time. Theoperation performing back-propagation is expressed by the adjoint formof the above given integral Equations (9)-(12) in the frequency domainor (15)-(18) in the time domain, which is obtained by swapping the signof iω or the time delay function τ(x, x^(R)), respectively. Thus, in thefirst embodiment of the method of the invention, the up-going anddown-going pressure fields at a point {tilde over (x)} below x^(R) aregiven by

$\begin{matrix}{{{p^{up}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\omega}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{p^{up}\left( {x^{R},\omega} \right)}{\exp \left( {{\omega\tau}\left( {\overset{\sim}{x},x^{R}} \right)} \right)}\ {S}}}}}{and}} & (19) \\{{p^{dn}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\omega}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{p^{dn}\left( {x^{R},\omega} \right)}{\exp \left( {{\omega\tau}\left( {\overset{\sim}{x},x^{R}} \right)} \right)}\ {S}}}}} & (20)\end{matrix}$

in the frequency domain and by

$\begin{matrix}{{{p^{up}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{p^{up}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}}{and}} & (21) \\{{p^{dn}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{p^{dn}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}} & (22)\end{matrix}$

in the time domain. The corresponding equations for up-going anddown-going normal components of the particle velocity field read:

$\begin{matrix}{{{v_{n}^{up}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\omega}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{v_{n}^{up}\left( {x^{R},\omega} \right)}{\exp \left( {{\omega\tau}\left( {\overset{\sim}{x},x^{R}} \right)} \right)}\ {S}}}}}{and}} & (23) \\{{v_{n}^{dn}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\omega}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{v_{n}^{dn}\left( {x^{R},\omega} \right)}{\exp \left( {{\omega\tau}\left( {\overset{\sim}{x},x^{R}} \right)} \right)}\ {S}}}}} & (24)\end{matrix}$

in the frequency domain and

$\begin{matrix}{{{v_{n}^{up}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{v_{n}^{up}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}}{and}} & (25) \\{{v_{n}^{dn}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\pi}}{\int_{S}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{v_{n}^{dn}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}} & (26)\end{matrix}$

in the time domain. For {tilde over (x)}=x in Equations (19) through(26), the output is given at the original acquisition surface S. Again,either the time-domain extrapolation expressions in Equations (21)-(22)and (25)-(26) or the frequency-domain extrapolation expressions inEquations (19)-(20) and (23)-(24) can be used for performing this seconddatuming step in this first embodiment of the invention.

FIGS. 3-5 and 7 show flowcharts illustrating embodiments of theinvention for wavefield separation for dual-sensor data. FIG. 3 is aflowchart illustrating a general embodiment of the invention. FIGS. 4-5and 7 are flowcharts further illustrating more particular embodiments ofthe invention as described in FIG. 3. FIGS. 1-2 and 6 are views andsections that illustrate aspects of the invention as described in theflowcharts in FIGS. 3-5 and 7.

FIG. 3 is a flowchart illustrating an embodiment of the invention forwavefield separation for dual-sensor data. At block 30, a first weightedintegral operator is applied to dual-sensor data to extrapolate thedual-sensor data to a first position above the acquisition surface,generating extrapolated data. The first weighted integral operator canbe applied in the time-space domain or the frequency-space domain.

At block 31, a second weighted integral operator, is applied to theextrapolated data from block 30 to extrapolate the extrapolated data toa second position, generating wavefield separated data, wherein one ofthe integral operators is applied to a scaled combination of the dualsensor data. The second weighted integral operator can be applied in thetime-space domain or the frequency-space domain.

FIG. 4 is a flowchart illustrating a first embodiment of the inventionfor wavefield separation for dual-sensor data.

At block 40, common-shot records of dual-sensor data are acquired at anacquisition surface below a sea surface.

At block 41, a first position level is selected above the acquisitionsurface from block 40.

At block 42, the dual sensor data are extrapolated from the acquisitionsurface from block 40 to the first position level from block 41, byapplying a first weighted integral operator to a scaled combination ofthe dual-sensor data from block 40, generating extrapolated up-going anddown-going data. The first weighted integral operator can be applied inthe time-space domain or the frequency-space domain.

At block 43, a second position level is selected. The second positionlevel can be the original acquisition surface from block 40, the seasurface from block 40, or any other position level that is useful toselect.

At block 44, the extrapolated up-going and down-going data from block 42are extrapolated from the first position level from block 41 to thesecond position level from block 43, by applying a second weightedintegral operator, generating wavefield separated data. The secondweighted integral operator can be applied in the time-space domain orthe frequency-space domain.

The second embodiment of the method of the invention is wavefieldseparation using vector-weighted Kirchhoff-type datuming. The secondembodiment uses only one recorded wavefield in the first datuming step,either the hydrophone record or the geophone record. In this datumingstep, the recorded wavefield is appropriately scaled. The same recordedwavefield is datumed without the scaling. Then, both records areback-propagated to the original acquisition surface. The ratio betweenthe scaled record and the unscaled record yields an estimate of theobliquity factor needed for wavefield separation. One of the hydrophoneor geophone records is scaled by the obliquity factor estimate and thenadded to or subtracted from (as appropriate) the other record to yieldup-going or down-going fields, respectively. The obliquity factorestimate is an amplitude weighted average of obliquity factorsassociated with several crossing events at a particular time sample andis given as a function of space and time. The output of this secondembodiment is, again, wavefield separated common-shot ensembles ofseismic traces.

In this second embodiment, only one of the recorded wavefields is usedin the first datuming step corresponding to Equations (9)-(12) or(15)-(18). In the following discussion, the method is described in termsof the pressure field for illustrative purposes only. The invention isnot restricted by this choice of field. The first datuming step isperformed with a two-component vector of additional weight functions andcreates a two-component vector of output pressure fields. The vector ofweight functions is given by:

w ₁(x)=(cos[θ(x,x ^(R))],1)^(T).   (27)

Using the vector of weight functions from Equation (27), the outputvector of pressure fields is given by:

$\begin{matrix}{{{p_{1}\left( {x^{R},t} \right)} \cong {\frac{1}{2\pi}{\int_{S}{{W\left( {x,x^{R}} \right)}{w_{1}(x)}\frac{\partial}{\partial t}{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}\ {S}}}}},} & (28)\end{matrix}$

where p₁(x^(R),t) designates the two-component vector:

p ₁(x ^(R) ,t)=(p ₁ ^(w)(x ^(R) ,t),p ₁(x ^(R) ,t)^(T).   (29)

Each of the two components of this vector of pressure fields p₁(x^(R),t)in Equation (29) is datumed back to the original acquisition surface,yielding a component of a second vector of pressure fields. These twosecond datuming equations are given by:

$\begin{matrix}{{{p_{2}^{w}\left( {x,t} \right)} \cong {{- \frac{1}{2\pi}}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{\partial}{\partial t}{p_{1}^{w}\left( {x^{R},{t + {\tau \left( {x,x^{R}} \right)}}} \right)}\ {S}}}}}{and}} & (30) \\{{p_{2}\left( {x,t} \right)} \cong {{- \frac{1}{2\pi}}{\int_{S}{{W\left( {x,x^{R}} \right)}\frac{\partial}{\partial t}{p_{1}\left( {x^{R},{t + {\tau \left( {x,x^{R}} \right)}}} \right)}\ {{S}.}}}}} & (31)\end{matrix}$

The ratio of the two pressure fields determined in Equations (30) and(31) yields the obliquity factor as a function of space and time:

$\begin{matrix}{{\cos\left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack} = {\frac{p_{2}^{w}\left( {x,t} \right)}{p_{2}\left( {x,t} \right)}.}} & (32)\end{matrix}$

For stability reasons and to avoid division by zero when p₂(x,t)=0 inthe denominator in Equation (32), it is advantageous to build theenvelope of the two pressure fields prior to the division. Thus, inanother embodiment of the invention, the obliquity factor is given by:

$\begin{matrix}{{\cos\left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack} = {\frac{{env}\left\lbrack {p_{2}^{w}\left( {x,t} \right)} \right\rbrack}{{env}\left\lbrack {p_{2}\left( {x,t} \right)} \right\rbrack}.}} & (33)\end{matrix}$

The term cos└{circumflex over (θ)}(x,t)┘ is used instead of cos[θ(x,t)]for the obliquity factor in Equations (32) and (33) to reflect that anamplitude weighted average of individual obliquity factors is consideredin the case of conflicting dips. This procedure provides the obliquityfactor that will be used in this second embodiment of the method of theinvention.

The case of three crossing events at a particular time sample andparticular position will illustrate this process of taking a weightedaverage of obliquity factors in more detail. The two pressure fieldsfrom Equations (30) and (31), obtained after the two datuming steps, arethen given by

p ₂ ^(w)(x ₀ ,t ₀)=a cos[θ_(a)(x ₀ ,t ₀)]+b cos[θ_(b)(x ₀ ,t ₀)]+ccos[θ_(c)(x ₀ ,t ₀)]  (34)

and

p ₂(x ₀ ,t ₀)=a+b +c,   (35)

where a, b, and c are the amplitudes of the three crossing events.

Consequently,

$\begin{matrix}{{{\cos \left\lbrack {\hat{\theta}\left( {x_{0},t_{0}} \right)} \right\rbrack} = \frac{\begin{matrix}{{a\; {\cos \left\lbrack {\theta_{a}\left( {x_{0},t_{0}} \right)} \right\rbrack}} + {b\; {\cos \left\lbrack {\theta_{b}\left( {x_{0},t_{0}} \right)} \right\rbrack}} +} \\{c\; {\cos \left\lbrack {\theta_{c}\left( {x_{0},t_{0}} \right)} \right\rbrack}}\end{matrix}}{a + b + c}},} & (36)\end{matrix}$

which is a weighted average of all three obliquity factors with theamplitudes of the individual events as weights.

Using the obliquity factor from either Equation (32) or (33), theseparated wavefields are then calculated using the following equations:

$\begin{matrix}{{{p^{up}\left( {x,t} \right)} = {\frac{1}{2}\left( {{p\left( {x,t} \right)} - {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{v_{n}\left( {x,t} \right)}}} \right)}},} & (37) \\{{{p^{dn}\left( {x,t} \right)} = {\frac{1}{2}\left( {{p\left( {x,t} \right)} + {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{v_{n}\left( {x,t} \right)}}} \right)}},} & (38) \\{{{v_{n}^{up}\left( {x,t} \right)} = {\frac{1}{2}\left( {{v_{n}\left( {x,t} \right)} - {\frac{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,t} \right)}}} \right)}},{and}} & (39) \\{{v_{n}^{dn}\left( {x,t} \right)} = {\frac{1}{2}{\left( {{v_{n}\left( {x,t} \right)} + {\frac{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,t} \right)}}} \right).}}} & (40)\end{matrix}$

FIG. 5 is a flowchart illustrating the second embodiment of theinvention for wavefield separation for dual-sensor data.

At block 50, common-shot records of dual-sensor data are acquired at anacquisition surface below a sea surface.

At block 51 a first position level above the acquisition surface fromblock 50 is selected.

At block 52 one selected field of the dual sensor data from block 50 isextrapolated from the acquisition surface from block 50 to the firstposition level from block 52, by applying a first weighted integraloperator to a scaled two-component vector of the one selected field,generating extrapolated two-component data. The first weighted integraloperator can be applied in the time-space domain or the frequency-spacedomain.

At block 53, the extrapolated two-component data from block 52 areextrapolated from the first position level from block 51 to theacquisition surface from block 50, by applying a second weightedintegral operator, generating two-component field data. The secondweighted integral operator can be applied in the time-space domain orthe frequency-space domain. However, if the results are in thefrequency-space domain, then the results must be transformed to thetime-space domain before the division in block 54.

At block 54, a first component is divided by a second component of thetwo-component field data from block 53, generating an estimatedobliquity factor.

At block 55, the estimated obliquity factor from block 54 is used as ascaler in scaled combinations of the dual-sensor data, generatingwavefield separated data.

The third embodiment of the method of the invention is essentially aKirchhoff-type depth migration with inherent wavefield separation. Thisthird embodiment is similar to the first embodiment in that the samecombination of common-shot records of the recorded wavefields is used.The recorded wavefields are not, however, extrapolated to a level abovethe acquisition surface, but migrated into the depth domain instead. Theoutput of this third embodiment is depth migrated volumes and/orcommon-image gathers (CIG).

The basic principle behind Kirchhoff-type depth migration is to sum allthe samples from traces within a given aperture around the lateraloutput position into the same output sample for which the recording timeequals the traveltimes along rays connecting the source positions ofthese traces with their receiver positions through the output locationunder consideration. The region where the traveltime operator is tangentto the seismic event in the pre-stack data contributes the major part tothe output amplitude. The principle of stationary phase thusautomatically selects the correct obliquity factor. The traveltimeoperator is usually called diffraction traveltime function as itdescribes the traveltimes of a seismic event related to a diffractionpoint in depth. Therefore, this migration is sometimes also calleddiffraction stack migration. The traveltime operator will be denoted byτ_(D) in the following.

Kirchhoff depth migration is expressed as a weighted summation along thediffraction traveltime operator:

$\begin{matrix}{{{I(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{\partial{p\left( {x^{S},x^{G},t} \right)}}{\partial t}}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},} & (41)\end{matrix}$

where M denotes the output depth location (the image location) underconsideration, I is the output amplitude (the image), A is the migrationaperture, and W (x^(S),x^(G),M) is a weighting factor accounting forgeometrical spreading and phase changes due to caustics. The inputpre-stack pressure data is given by p(x^(S),x^(G),t). An analogousexpression for the normal component of the particle velocity fieldυ_(n)(x^(S),x^(G),t) is obtained by replacing p(x^(S),x^(G),t) byυ_(n)(x^(S),x^(G),t) in Equation (41). The integral above is formulatedusing shot (x^(S)) and receiver (x^(G)) coordinates, for illustrativepurposes only. The choice of coordinates is not intended to be alimitation of the invention. There exist alternative formulations indifferent coordinate systems and everything following from now onapplies equally well to these alternative formulations.

Analogous to the datuming in the previous section, the Kirchhoffmigration in the above stated form is based on a high-frequencyapproximation. Thus, the input data can be treated as being composed ofseismic events related to locally plane waves. Keeping the sourceposition fixed and taking the first derivative with respect to thereceiver position of the diffraction traveltime function provides thetime dip of the traveltime operator. Using the local velocity at thereceiver position, this time dip can be used to calculate the equivalentto cos└θ(x,x^(R))┘ appearing in the time domain datuming approach, i.e.the cosine of the emergence angle of the ray connecting the imagelocation M with the receiver location x^(G). This is denoted bycos└θ(x^(G),M)┘ in the following. This angle is more directly obtainedfrom ray-tracing, a common method to calculate the diffractiontraveltime operator. Using the following abbreviations

$\begin{matrix}{{{p_{1}^{G}\left( {x^{S},x^{G},M} \right)} = \frac{\partial{\tau_{D}\left( {x^{S},x^{G},M} \right)}}{\partial x_{1}^{G}}}{and}} & (42) \\{{{p_{2}^{G}\left( {x^{S},x^{G},M} \right)} = \frac{\partial{\tau_{D}\left( {x^{S},x^{G},M} \right)}}{\partial x_{2}^{G}}},{\cos \left\lfloor {\theta \left( {x^{G},M} \right)} \right\rfloor \mspace{14mu} {is}\mspace{14mu} {obtained}\mspace{14mu} {by}}} & (43) \\{{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack} = {{c\left( x^{G} \right)}\sqrt{\frac{1}{c^{2}\left( x^{G} \right)} - \left( p_{1}^{G} \right)^{2} - \left( p_{2}^{G} \right)^{2}}}} & (44)\end{matrix}$

The geometrical situation just described is depicted in FIG. 6, which isa schematic side view of a single selected ray corresponding to a timesample from a selected trace which is to be migrated. Analogous to thetime domain datuming procedure it is possible to perform the wavefieldseparation during Kirchhoff depth migration. This migration isaccomplished by the following integral equations

$\begin{matrix}{{{I_{p}^{up}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{p\left( {x^{S},x^{G},t} \right)} - {\frac{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{v_{n}\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},} & (45) \\{{{I_{p}^{dn}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{p\left( {x^{S},x^{G},t} \right)} + {\frac{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{v_{n}\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},} & (46) \\{{{I_{v}^{up}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{v_{n}\left( {x^{S},x^{G},t} \right)} - {\frac{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{p\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},\mspace{79mu} {and}} & (47) \\{{I_{v}^{dn}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2}\frac{\partial}{\partial t}\left( {{v_{n}\left( {x^{S},x^{G},t} \right)} + {\frac{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{p\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}{{A}.} \right.} & (48)\end{matrix}$

The principle of stationary phase automatically selects the correctscaling of the particle velocity field or the pressure field inEquations (45)-(48), as only contributions from a small zone around thestationary point constructively contribute to the integration.

FIG. 7 is a flowchart illustrating an embodiment of the invention forwavefield separation for dual-sensor data.

At block 70, common-shot records of dual-sensor data are acquired at anacquisition surface.

At block 71, the dual-sensor data from block 70 are migrated from theacquisition surface from block 70 into a depth domain, by applying aweighted integral operator to a scaled combination of the dual-sensordata from block 70 in the time-space domain, generating wavefieldseparated data.

For a practical implementation, all integrations in the embodimentsdescribed above have to be replaced by discrete summations. Areaelements, which appear inside the integral expressions, are given by amultiplication of the receiver spacings in both lateral directions for aregular grid of receivers. In any case, including irregular geometries,these surface elements can be calculated by means well-known to peopleskilled in the art. The same applies to the weights appearing in theintegral expressions. They are also well-known in the art. Theimplementations can be based on standard implementations ofKirchhoff-type datuming and migration schemes.

It should be understood that the preceding is merely a detaileddescription of specific embodiments of this invention and that numerouschanges, modifications, and alternatives to the disclosed embodimentscan be made in accordance with the disclosure here without departingfrom the scope of the invention. The preceding description, therefore,is not meant to limit the scope of the invention. Rather, the scope ofthe invention is to be determined only by the appended claims and theirequivalents.

1. A method for wavefield separation for dual-sensor data, comprising:applying a first weighted integral operator to dual-sensor data toextrapolate the dual-sensor data to a first position above anacquisition surface, generating extrapolated data; and applying a secondweighted integral operator to the extrapolated data to extrapolate theextrapolated data to a second position, generating wavefield separateddata, wherein one of the integral operators is applied to a scaledcombination of the dual sensor data.
 2. The method of claim 1, whereinthe applying a first weighted integral operator comprises: acquiringcommon-shot records of dual-sensor data at an acquisition surface belowa sea surface; selecting a first position level above the acquisitionsurface; and extrapolating the dual sensor data in the time-domain fromthe acquisition surface to the first position level, by applying a firstweighted integral operator to a scaled combination of the dual-sensordata, generating extrapolated up-going and down-going data.
 3. Themethod of claim 2, wherein the extrapolating the dual-sensor data isaccomplished by applying the following equations in the frequency-spacedomain:${{p^{up}\left( {x^{R},\omega} \right)} \cong {\frac{\; \omega}{2\; \pi} {\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \left( {{p\left( {x,\omega} \right)} - {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,\omega} \right)}}} \right){\exp \left\lbrack {{- }\; \omega \; {\tau \left( {x,x^{R}} \right)}} \right\rbrack}\ {S}}}}},{{p^{dn}\left( {x^{R},\omega} \right)} \cong {\frac{\; \omega}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\left( {{p\left( {x,\omega} \right)} + {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,\omega} \right)}}} \right){\exp \left\lbrack {{- }\; \omega \; {\tau \left( {x,x^{R}} \right)}} \right\rbrack}\ {S}}}}},{{v_{n}^{up}\left( {x^{R},\omega} \right)} \cong {\frac{\; \omega}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\left( {{v_{n}\left( {x,\omega} \right)} - {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,\omega} \right)}}} \right){\exp \left( {{- }\; \omega \; {\tau \left( {x,x^{R}} \right)}} \right)}\ {S}}}}},\mspace{79mu} {{{and}{v_{n}^{dn}\left( {x^{R},\omega} \right)}} \cong {\frac{\; \omega}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2}\left( {{v_{n}\left( {x,\omega} \right)} + {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,\omega} \right)}}} \right){\exp \left( {{- }\; \omega \; {\tau \left( {x,x^{R}} \right)}} \right)}\ {S}}}}},$where p^(up) (x^(R),ω) is an up-going component of pressure extrapolatedto x^(R), p^(dn)(x^(R),ω) is a down-going component of pressureextrapolated to x^(R), υ_(n) ^(dn)(x^(R),ω) is an up-going component ofnormal component of particle velocity extrapolated to x^(R), υ_(n)^(dn)(x^(R),ω) is a down-going component of normal component of particlevelocity extrapolated to x^(R), x is a point on an acquisition surfaceS, x^(R) is a point above the acquisition surface S, ω is angularfrequency, W (x,x^(R)) is a weighting factor, p(x,ω) is pressuremeasured at x, υ_(n)(x,ω) is a normal component of particle velocitymeasured at x, ρ(x) is density, c(x) is speed of sound, θ(x,x^(R)) is anangle between a ray connecting x and x^(R) and a unit normal vector tothe acquisition surface S, and τ(x, x^(R)) is traveltime along the rayconnecting x and x^(R).
 4. The method of claim 2, wherein theextrapolating the dual-sensor data is accomplished by applying thefollowing equations in the time-space domain:${{p^{up}\left( {x^{R},t} \right)} \cong {\frac{1}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \frac{\partial}{\partial t} \left( {{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)} - {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right)\ {S}}}}},{{p^{dn}\left( {x^{R},t} \right)} \cong {\frac{1}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \frac{\partial}{\partial t} \left( {{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)} + {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right)\ {S}}}}},{{v_{n}^{up}\left( {x^{R},t} \right)} \cong {\frac{1}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \frac{\partial}{\partial t} \left( {{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)} - {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right)\ {S}}}}},\mspace{79mu} {{{and}{v_{n}^{dn}\left( {x^{R},t} \right)}} \cong {\frac{1}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{1}{2} \frac{\partial}{\partial t} \left( {{v_{n}\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)} + {\frac{\cos \left\lbrack {\theta \left( {x,x^{R}} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}}} \right)\ {S}}}}},$where p^(up)(x^(R),t) is an up-going component of pressure extrapolatedto x^(R), p^(dn)(x^(R),t) is a down-going component of pressureextrapolated to x^(R), υ_(n) ^(up)(x^(R),t) is an up-going component ofnormal component of particle velocity extrapolated to x^(R), υ_(n)^(dn)(x^(R),t) is a down-going component of normal component of particlevelocity extrapolated to x^(R), x is a point on an acquisition surfaceS, x^(R) is a point above the acquisition surface S, t is time,W(x,x^(R) )is a weighting factor, p(x,t) is pressure measured at x,v_(n)(x,t) is a normal component of particle velocity measured at x,ρ(x) is density, c(x) is speed of sound, θ(x,x^(R)) is an angle betweena ray connecting x and x^(R) and a unit normal vector to the acquisitionsurface S, and τ(x,x^(R)) is traveltime along the ray connecting x andx^(R).
 5. The method of claim 2, wherein the applying a second weightedintegral operator comprises: selecting a second position level; andextrapolating the extrapolated up-going and down-going data from thefirst position level to the second position level, by applying a secondweighted integral operator generating wavefield separated data.
 6. Themethod of claim 5, wherein the extrapolating the extrapolated up-goingand down-going data is accomplished by applying the following equationsin the frequency-space domain:${{p^{up}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\; \omega}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{p^{up}\left( {x^{R},\omega} \right)}{\exp \left( {\; \omega \; {\tau \left( {\overset{\sim}{x},x^{R}} \right)}} \right)}\ {S}}}}},{{p^{dn}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\; \omega}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{p^{dn}\left( {x^{R},\omega} \right)}{\exp \left( {\; \omega \; {\tau \left( {\overset{\sim}{x},x^{R}} \right)}} \right)}\ {S}}}}},{{v_{n}^{up}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\; \omega}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{v_{n}^{up}\left( {x^{R},\omega} \right)}{\exp \left( {\; \omega \; {\tau \left( {\overset{\sim}{x},x^{R}} \right)}} \right)}\ {S}}}}},{and}$${{v_{n}^{dn}\left( {\overset{\sim}{x},\omega} \right)} \cong {{- \frac{\; \omega}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}{v_{n}^{dn}\left( {x^{R},\omega} \right)}{\exp \left( {\; \omega \; {\tau \left( {\overset{\sim}{x},x^{R}} \right)}} \right)}\ {S}}}}},$where p^(up)({tilde over (x)},ω) is an up-going component of pressureextrapolated to {tilde over (x)}, p^(dn)({tilde over (x)},ω) is adown-going component of pressure extrapolated to {tilde over (x)}, ν_(n)^(up)({tilde over (x)},ω) is an up-going component of normal componentof particle velocity extrapolated to {tilde over (x)}, ν_(n)^(dn)({tilde over (x)},ω) is a down-going component of normal componentof particle velocity extrapolated to {tilde over (x)}, x^(R) is a pointabove an acquisition surface S, {tilde over (x)} is a point below x^(R),ω is angular frequency, W(x^(R),{tilde over (x)}) is a weighting factor,p^(up)(x^(R),ω)) is an up-going component of pressure extrapolated tox^(R), ν_(n) ^(up)(x^(R),ω)) is a down-going component of pressureextrapolated to x^(R), ν_(n) ^(up)(x^(R),ω) is an up-going component ofnormal component of particle velocity extrapolated to x^(R), v_(n)^(dn)(x^(R),ω) is a down-going component of normal component of particlevelocity extrapolated to x^(R), and τ({tilde over (x)},x^(R)) istraveltime along the ray connecting {tilde over (x)} and x^(R).
 7. Themethod of claim 5, wherein the extrapolating the extrapolated up-goingand down-going data is accomplished by applying the following equationsin the time-space domain:${{p^{up}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{p^{up}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}},{{p^{dn}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{p^{dn}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}},{{v_{n}^{up}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{v_{n}^{up}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}},{and}$${{v_{n}^{dn}\left( {\overset{\sim}{x},t} \right)} \cong {{- \frac{1}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x^{R},\overset{\sim}{x}} \right)}\frac{\partial}{\partial t}{v_{n}^{dn}\left( {x^{R},{t + {\tau \left( {\overset{\sim}{x},x^{R}} \right)}}} \right)}\ {S}}}}},$where p^(up)({tilde over (x)},t) is an up-going component of pressureextrapolated to {tilde over (x)}, p^(dn)({tilde over (x)},t) is adown-going component of pressure extrapolated to {tilde over (x)}, υ_(n)^(up)({tilde over (x)},t) is an up-going component of normal componentof particle velocity extrapolated to {tilde over (x)}, υ_(n)^(dn)({tilde over (x)},t) is a down-going component of normal componentof particle velocity extrapolated to {tilde over (x)}, x^(R) is a pointabove an acquisition surface S, {tilde over (x)} is a point below x^(R),t is time, W(x^(R),{tilde over (x)}) is a weighting factor,p^(up)(x^(R),t) is an up-going component of pressure extrapolated tox^(R), p^(dn)(x^(R),t) is a down-going component of pressureextrapolated to x^(R), υ_(n) ^(up)(x^(R),t) is an up-going component ofnormal component of particle velocity extrapolated to x^(R), υ_(n)^(dn)(x^(R),t) is a down-going component of normal component of particlevelocity extrapolated to x^(R), and τ({tilde over (x)}, x^(R)) istraveltime along the ray connecting {tilde over (x)} and x^(R).
 8. Themethod of claim 1, wherein the applying a first weighted integraloperator comprises: acquiring common-shot records of dual-sensorstreamer data at an acquisition surface below a sea surface; selecting afirst position level above the acquisition surface; and extrapolatingone selected field of the dual sensor data in the time-domain from theacquisition surface to the first position level, by applying a firstweighted integral operator to a scaled two-component vector of the oneselected field, generating extrapolated two-component data.
 9. Themethod of claim 8, wherein the extrapolating one selected field of thedual sensor data is accomplished by applying the following equations inthe time-space domain:${{p_{1}\left( {x^{R},t} \right)} \cong {\frac{1}{2\; \pi}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}{w_{1}(x)}\frac{\partial}{\partial t}{p\left( {x,{t - {\tau \left( {x,x^{R}} \right)}}} \right)}\ {S}}}}},$where p₁(x^(R),t) designates a two-component pressure vector:p ₁(x ^(R) ,t)=(p ₁ ^(w)(x ^(R) ,t),p ₁(x ^(R) ,t)^(T), W (x, x^(R)) isa weighting factor, w₁(x) is a vector of weight functions given by:w ₁(x)=(cos[θ(x,x ^(R))],1)^(T), x is a point on an acquisition surfaceS, x^(R) is a point above the acquisition surface S, p(x,t) is pressuremeasured at x on the acquisition surface S, p_(x) ^(w)(x^(R),t) is afirst weighted pressure component at the point x^(R), p₁(x^(R),t) is afirst unweighted pressure component at the point x^(R), θ(x,x^(R)) is anangle between a ray connecting x and x^(R) and a unit normal vector tothe acquisition surface S, and τ(x,x^(R)) is traveltime along a rayconnecting x and x^(R).
 10. The method of claim 8, wherein the applyinga second weighted integral operator comprises: extrapolating theextrapolated two-component data from the first position level to theacquisition surface, by applying a second weighted integral operator,generating two-component field data; dividing a first component by asecond component of the two-component field data, generating anestimated obliquity factor; and using the estimated obliquity factor asa scaler in scaled combinations of the dual-sensor data, generatingwavefield separated data.
 11. The method of claim 10, wherein theextrapolating the extrapolated two-component data is accomplished byapplying the following equations in the time-space domain:${p_{2}^{w}\left( {x,t} \right)} \cong {{- \frac{1}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{\partial}{\partial t}{p_{1}^{w}\left( {x^{R},{t + {\tau \left( {x,x^{R}} \right)}}} \right)}\ {S}}}}$and${p_{2}\left( {x,t} \right)} \cong {{- \frac{1}{2\; \pi}}{\int_{S}^{\;}{{W\left( {x,x^{R}} \right)}\frac{\partial}{\partial t}{p_{1}\left( {x^{R},{t + {\tau \left( {x,x^{R}} \right)}}} \right)}\ {{S}.}}}}$where p₂ ^(w)(x,t) is a second weighted pressure component at the pointx, p₂(x,t) is a second unweighted pressure component at the point x, W(x,x^(R)) is a weighting factor, x is a point on an acquisition surfaceS, x^(R) is a point above the acquisition surface S, p₁ ^(w)(x^(R),t) isa first weighted pressure component at the point x^(R), p₁(x^(R),t) is afirst unweighted pressure component at the point x^(R), and τ(x,x^(R))is traveltime along a ray connecting x and x^(R).
 12. The method ofclaim 11, wherein the dividing a first component by a second componentis accomplished by applying the following equation in the time-spacedomain:${\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack} = {\frac{p_{2}^{w}\left( {x,t} \right)}{p_{2}\left( {x,t} \right)}.}$where cos└{tilde over (θ)}(x,t)┘ is an estimated obliquity factor. 13.The method of claim 11, wherein the dividing a first component by asecond component is accomplished by applying the following equation inthe time-space domain:${\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack} = {\frac{{env}\left\lbrack {p_{2}^{w}\left( {x,t} \right)} \right\rbrack}{{env}\left\lbrack {p_{2}\left( {x,t} \right)} \right\rbrack}.}$where cos└{tilde over (θ)}(x,t)┘ is an estimated obliquity factor. 14.The method of claim 10, wherein the using the estimated obliquity factoras a scaler is accomplished by applying the following equations in thetime-space domain:${{p^{up}\left( {x,t} \right)} = {\frac{1}{2}\left( {{p\left( {x,t} \right)} - {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{v_{n}\left( {x,t} \right)}}} \right)}},{{p^{dn}\left( {x,t} \right)} = {\frac{1}{2}\left( {{p\left( {x,t} \right)} + {\frac{{\rho (x)}{c(x)}}{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{v_{n}\left( {x,t} \right)}}} \right)}},{{v_{n}^{up}\left( {x,t} \right)} = {\frac{1}{2}\left( {{v_{n}\left( {x,t} \right)} - {\frac{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,t} \right)}}} \right)}},{and}$${v_{n}^{dn}\left( {x,t} \right)} = {\frac{1}{2}{\left( {{v_{n}\left( {x,t} \right)} + {\frac{\cos \left\lbrack {\hat{\theta}\left( {x,t} \right)} \right\rbrack}{{\rho (x)}{c(x)}}{p\left( {x,t} \right)}}} \right).}}$where p^(up)(x,t) is an up-going component of pressure extrapolated tox, p^(dn)(x,t) is a down-going component of pressure extrapolated to x,υ_(n) ^(up)(x,t) is an up-going component of normal component ofparticle velocity extrapolated to x, υ_(n) ^(dn)(x,t) is a down-goingcomponent of normal component of particle velocity extrapolated to x,x^(R) is a point above an acquisition surface S, p(x,t) is pressuremeasured at x, υ_(n)(x,t) is a normal component of particle velocitymeasured at x, ρ(x) is density, c(x) is speed of sound, and cos└{tildeover (θ)}(x,t)┘ is an estimated obliquity factor.
 15. The method ofclaim 1, wherein the applying a first weighted integral operatorcomprises: acquiring common-shot records of dual-sensor data at anacquisition surface.
 16. The method of claim 15, wherein the applying asecond weighted integral operator comprises: migrating the dual-sensordata from the acquisition surface into a depth domain, by applying inthe time-space domain a weighted integral operator to a scaledcombination of the dual-sensor data, generating wavefield separateddata.
 17. The method of claim 16, wherein the migrating the dual-sensordata is accomplished by applying the following equations:${{I_{p}^{up}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2} \frac{\partial}{\partial t} \left( {{p\left( {x^{S},x^{G},t} \right)} - {\frac{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{v_{n}\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},{{I_{p}^{dn}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2}\frac{\partial}{\partial t} \left( {{p\left( {x^{S},x^{G},t} \right)} + {\frac{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{v_{n}\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},{{I_{v}^{up}(M)} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2}\frac{\partial}{\partial t} \left( {{v_{n}\left( {x^{S},x^{G},t} \right)} - {\frac{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{p\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}\ {A} \right.},\mspace{79mu} {{{and}{I_{v}^{dn}(M)}} = \left. {{- \frac{1}{2\; \pi}}{\int_{A}^{\;}{{W\left( {x^{S},x^{G},M} \right)}\frac{1}{2} \frac{\partial}{\partial t} \left( {{v_{n}\left( {x^{S},x^{G},t} \right)} + {\frac{\cos \left\lbrack {\theta \left( {x^{G},M} \right)} \right\rbrack}{{\rho \left( x^{G} \right)}{c\left( x^{G} \right)}}{p\left( {x^{S},x^{G},t} \right)}}} \right)}}} \middle| {}_{t = {\tau_{D}{({x^{S},x^{G},M})}}}{A} \right.},$where I_(p) ^(up)(M) is an up-going amplitude of pressure at imagelocation M, I_(p) ^(dn)(M) is a downgoing amplitude of pressure at imagelocation M, I_(υ) ^(up)(M) is an up-going amplitude of normal particlevelocity at image location M, I_(υ) ^(dn)(M) is a down-going amplitudeof normal particle velocity at image location M, M is an image location,A is a migration aperture, W(x^(S),x^(G),M) is a weighting factoraccounting for geometrical spreading and phase changes due to caustics,x^(S) is a source location, x^(G) is a receiver location, t is time,p(x^(S),x^(G),t) is pressure, υ_(n)(x^(S),x^(G),t) is a normal componentof the particle velocity field, ρ(x^(G)) is density, c(x^(G)) is speedof sound, θ(x^(G),M) is emergence angle of the ray connecting the imagelocation M with the receiver location x^(G), and τ_(D)(x^(S),x^(G),M) isdiffraction traveltime operator.